/* 
 * File:   SVD.h
 * Author: gustavo
 *
 * Created on August 22, 2012, 10:40 AM
 */

#ifndef SVD_HPP
#define	SVD_HPP

#include "../matrix.hpp"

#include "../impl/mdec/svd.hpp"

namespace imp {

	namespace math {

		template <typename T> class SVD {
			////////

		public:
			Matrix<T> U;
			Matrix<T> S;
			Matrix<T> V;

		public: /* Constructors */
			SVD(Matrix<T>& orig);
			virtual ~SVD();

		public: /* Public Methods */
			const Matrix<T>& getU() const;
			const Matrix<T>& getS() const;
			const Matrix<T>& getV() const;

		};

		template <typename T> SVD<T>::SVD(Matrix<T>& M) {

			if (M.rows >= M.columns) {

				U.setSize(M.rows, M.columns);
				S.setSize(M.columns, 1);
				V.setSize(M.columns, M.columns);

				impl::svd(M.rows, M.columns, 1, 1, M.data, S.data, U.data, V.data);
			} else {
				
				Matrix<T> Mt;
				
				Mt = M.transpose();
				
				V.setSize(Mt.rows, Mt.columns);
				S.setSize(Mt.columns, 1);
				U.setSize(Mt.columns, Mt.columns);

				impl::svd(Mt.rows, Mt.columns, 1, 1, Mt.data, S.data, V.data, U.data);
			}

		}

		template <typename T> SVD<T>::~SVD() {

		}

		template <typename T> const Matrix<T>& SVD<T>::getS() const {
			return this->S;
		}

		template <typename T> const Matrix<T>& SVD<T>::getU() const {
			return this->U;
		}

		template <typename T> const Matrix<T>& SVD<T>::getV() const {
			return this->V;
		}


	}
}

#endif	/* SVD_HPP */

